Two-layer analytical model for estimation of layer thickness and flow using Diffuse Correlation Spectroscopy

Diffuse correlation spectroscopy (DCS) has been widely explored for its ability to measure cerebral blood flow (CBF), however, mostly under the assumption that the human head is homogenous. In addition to CBF, knowledge of extracerebral layers, such as skull thickness, can be informative and crucial for patient with brain complications such as traumatic brain injuries. To bridge the gap, this study explored the feasibility of simultaneously extracting skull thickness and flow in the cortex layer using DCS. We validated a two-layer analytical model that assumed the skull as top layer with a finite thickness and the brain cortex as bottom layer with semi-infinite geometry. The model fitted for thickness of the top layer and flow of the bottom layer, while assumed other parameters as constant. The accuracy of the two-layer model was tested against the conventional single-layer model using measurements from custom made two-layer phantoms mimicking skull and brain. We found that the fitted top layer thickness at each source detector (SD) distance is correlated with the expected thickness. For the fitted bottom layer flow, the two-layer model fits relatively consistent flow across all top layer thicknesses. In comparison, the conventional one-layer model increasingly underestimates the bottom layer flow as top layer thickness increases. The overall accuracy of estimating first layer thickness and flow depends on the SD distance in relationship to first layer thickness. Lastly, we quantified the influence of uncertainties in the optical properties of each layer. We found that uncertainties in the optical properties only mildly influence the fitted thickness and flow. In this work we demonstrate the feasibility of simultaneously extracting of layer thickness and flow using a two-layer DCS model. Findings from this work may introduce a robust and cost-effective approach towards simultaneous bedside assessment of skull thickness and cerebral blood flow.


Introduction
Microvascular blood flow ensures delivery of oxygen and nutrients to the tissue and subsequent removal of metabolic by-products from the tissue, hence crucial for healthy functionality disadvantages including poor safety profile as mentioned above [8][9][10]29]. In summary, continuous bedside monitoring to quantify both skull thickness and CBF is presently not available for clinical use. Towards this goal, we adapted a two-layer analytical model from the multi-layer model described by Li et al. [22]. In this model, we assumed the skull as top layer with finite thickness and the brain cortex as bottom layer with semi-infinite geometry; we then use it to fit for thickness in the top layer thickness and flow in the bottom layer. The flow extracted from the proposed model was tested against the conventional single-layer model using measurements from custom made two-layer phantoms mimicking skull and the brain for various source-detector (SD) separations. Lastly, since layer optical properties (absorption coefficient (μ a ) and reduced scattering coefficient (m 0 s )) of the phantoms were assumed in the model, an error analysis was conducted to explore how much an incorrect assumption of layer optical properties might alter the model performance. This work explores the feasibility of simultaneous extraction of layer thickness and flow, which may introduce a robust and cost-effective approach towards simultaneous bedside assessment of skull thickness and CBF.

Methods
In this section, the theory of the single-layer and two-layer analytical models for CW-DCS are explained. In addition, details on device instrumentation, data acquisition and processing strategies are described. Lastly, details on fabrication of the two-layer phantoms and measurement are outlined.

Single-layer analytical model
In DCS, the near-infrared (NIR) light is directed into a diffusive medium. Photons then undergo dynamic phase shifts from the moving scatterers, i.e., red blood cells (RBC) [2]. Such phase shifts give rise to temporally varying speckle patterns which contains dynamic information of the scatterers. DCS measures the temporal fluctuation of the reflected light, from which a metric for CBF can be computed. The temporal electric field autocorrelation function (G 1 ) of the scattered electric field (E(t)), which carries the dynamic properties of the scatterer, can be written as [30]: where the brackets hi denote the average over time t, and τ is the delay time.
In semi-infinite homogeneous media, G 1 can be described by the correlation diffusion equation: ρ is the source-detector distance, 1À R eff ; R eff ¼ À 1:44n À 2 þ 0:71n À 1 þ 0:064n þ 0:668, n is the refractive index, α describes the fraction of photon scattering events from moving scatterers such as red blood cells in the tissue, k 0 is the wave-number of the light in tissue, and hΔr 2 (τ)i is the mean square displacement of the moving scatterers. In most DCS experiments in living tissues, the hΔr 2 (τ)i has been found to be reasonably well approximated as an effective Brownian motion, i.e., hΔr 2 (τ)i = 6D B τ, where D B is the Brownian diffusion coefficient of the moving scatterers and thus the flow dependent parameter in the correlation-diffusion equation [30].
In practice, DCS measures the fluctuations of the scattered light intensity. The normalized temporal intensity autocorrelation function is calculated as where I is the measured intensity. The normalized temporal electric field autocorrelation function g 1 is related to g 2 by the Siegert relationship [2,11].
Another single-layer analytical model described by Li et al. will be introduced in the next session.

Multi-layer analytical model
Li et al. proposed a temporal field autocorrelation function with respect to a multi-layered tissue geometry [22]. In this manuscript, we adapted the model described by Li et al. and formulated it for a two-layer geometry to simulate the skull and the brain, as shown in Fig 1. The optical properties for each layer are represented by the transport mean free path l � n ¼ 1=m 0 s ðnÞ [31] and absorption mean free path l ðaÞ n ¼ 1=m ðnÞ a , where m ðnÞ a and m 0 s ðnÞ are absorption and reduced scattering coefficients for each layer with n = 1-2. The d n represents thicknesses for respective layers.
The temporal field autocorrelation function Gðr; tÞ ¼ hEðtÞ � E � ðt þ tÞi can be described by the correlation diffusion equation [12] as where s 0 is a point-like monochromatic light source located at r 0 = {ρ 0 = 0, z 0 } inside the first layer; ρ represents the transverse coordinate, and z 0 ¼ 1=m 0 s ð1Þ. The decorrelation due to the moving scatterers, α n , can be written as a n ðtÞ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 3 where for each layer, t ð0Þ n ¼ ðk 2 n D Bn Þ À 1 is the correlation time for a single scattering event. Here, k n is the wavenumber of the light at wavelength λ and D Bn is the Brownian diffusion coefficient for the n-th layer, a parameter that describes the scatterer dynamics inside each layer and thus can quantify the motion of RBC in thick/deep tissue [2,11]. The field autocorrelation at the surface, G 0 (r, τ), can be obtained by solving Eq (6) in the Fourier domain with respect to the transverse coordinate ρ asĜ where q is the radial spatial frequency. By applying the boundary conditions as described in Li et al., the Fourier transform of the field autocorrelation function of diffuse reflected light, G 0 (r, τ), measured at the surface can be written aŝ In this work, the human head was modeled as a two-layer medium, as shown in Fig 1. The first layer accounts for skull as well as the static (i.e., no moving scatterers) scalp and the second layer has a semi-infinite geometry that accounts for the brain cortex. For n = 2, d 2 ! 1: hence, the numerator and denominator forĜ 0 ðq; z ¼ 0; tÞ can be written as where z 0 0 ¼ 1=m 0 s ð1Þ, and for each layer, b n ðq; tÞ ¼ ½a 2 n ðtÞ þ q 2 � 1=2 ; D n ¼ cl � n =3 is the photon diffusion coefficient, and c is the speed of light. s are the absorption and reduced scattering coefficients, l � n is the transport mean free path, and l a n is the absorption mean free path, for the respective layers where n = 1-2. The top layer is the solid layer with finite thickness (d 1 ), the bottom layer is the liquid layer with semi-infinite geometry (d 2 ). https://doi.org/10.1371/journal.pone.0274258.g001 For n = 1, a homogenous layer, the numerator and denominator for Eq (9) are Lastly, by performing the inverse Fourier transform of Eq (9) with respect to q, the autocorrelation function for the electric field measured at the surface can be obtained as where J 0 is the zero-order Bessel function of the first kind.

DCS instrumentation and data analysis
A custom built DCS system was used in this work. Following that, an 8-channel hardware correlator (Flex05-8ch, Correlator.com, NJ), with 4 channels being available for each SD distance, was used to auto-correlate the detected intensities with the sampling frequency of 2 Hz, from which the intensity autocorrelation (g 2 ) was obtained. Time span for a single measurement was 100 s. Each phantom measurement was repeated ten times, between which the probe was taken off the phantom and put back in place. The photon count rate wasn't adjusted by tuning the laser power during the measurement. During data processing, the autocorrelation correction factor, β, was estimated from the first few data points (τ = 0.5×10 −6 s to τ = 0.8×10 −6 s) of the average of measured g 2 over time. The averaged g 2 was then converted to g 1 using Eq (5). Additionally, the SNR for each delay time over the 100 s measurement was calculated using SNR = (g 2 (τ)−1)/σ(g 2 (τ)) [17]. For the two-layer model, the g 1 curve from τ = 1×10 −6 s to τ = 3×10 −3 s for each measurement was then fitted to the solution of the correlation diffusion equation using Eq (14) with n = 2 to extract top layer thickness (d 1 ) and bottom layer Brownian diffusion coefficient (D B2 ) [2]. In addition, ground-truth D B2 was extracted using Eq (14) with n = 1 on the bottom-layer-only (homogenous) measurements. For comparison purpose, the g 1 curve above 0.7 was fitted to the conventional single-layer model (Eq (2)) to extract the Brownian diffusion coefficient times α (αD B ) [32].
For both one-parameter single-layer models-Eq (2) and Eq (14) with n = 1 -we used MATLAB (Mathworks, Inc., Natick, MA) "fminsearch" function to perform least-square fitting. For the two-parameter two-layer model (Eq (14) with n = 2), MATLAB "fmincon" function was used with SNR as a weighting function for each delay time. The weighting was implemented to reduce noise on the fitted curve. The lower and upper bound for d 1 and D B2 were set to be in the range of 0 to 20 mm and 0 to 10 −6 cm 2 /s respectively with the initialization set to 0 for both parameters. The D B1 was set to zero for the use of a solid and thus a static layer [22]. The photon diffusion coefficients, D 1 and D 2 , were calculated using measured optical properties and were held constant. These assumptions reduced the free parameters only to d 1 and D B2 during fitting. The average and standard deviation from the 10 repeated measurements were calculated for fitted d 1 and fitted D B2 .
To confirm the fitting quality of the two-layer model, we computed the residual sum of squares (RSS = ∑[g 1,fitted (τ)−g 1,measured (τ)] 2 ) between the fitted and measured g 1 and generated the corresponding contour plots with various d 1 and D B2 values. In addition, if the global minimum of RSS(d 1 , D B2 ) indicates a different d 1 and D B2 pair than the fitted values by MATLAB "fmincon", the d 1 and D B2 pair from the contour plot would be used as the fitted values, because "fmincon" does not guarantee convergence to the global minimum. To test how the uncertainty in optical properties will affect the two-layer model (Eq (14) with n = 2), we modified each of the μ a and m 0 s of each layer used in the fitting by ± 20%, and then performed the same two-parameter fitting procedure as described previously. Percent changes of the fitted d 1 from true d 1 were computed by 100%�(d 1,fitted −d 1,true )/d 1,true ; percent changes of the fitted D B2 from true D B2 (fitted by the true μ a and m 0 s using the Eq (14) with n = 1) were computed by 100%�(D B2,fitted −D B2,true )/D B2,true .

Two-layer phantom experiment
The top layer of the two-layer optical phantoms, as shown in Fig 2,  were added to the silicone base, mixed with a handheld electrical mixer, and then poured onto a petri dish. This recipe follows the standard mixing ratios as reported in the literature [33,34]. The thickness of the mixture inside the petri dish was varied to mimic different skull thicknesses. Following that, the mixture was degassed for~1 hour in a vacuum chamber to remove air bubbles and left to solidify for~24 hours. Lastly, the solid phantom was peeled out of the petri dish and the thickness (referred to as true thickness in later sections) was measured using a slide caliper. Human skull thickness varies with sex, age, etc., but it was reported to be within 1-10 mm [35,36]. As a result, in this work, seven different thickness values ranging from 2.01 mm to 8.08 mm were tested (P1 = 2.01, P2 = 3.28, P3 = 3.87, P4 = 4.65, P5 = 4.99, P6 = 5.96, and P7 = 8.08 mm) for the top layer.
The bottom layer in the two-layer phantoms was a liquid phantom which represented the human brain cortex, as shown in Fig 2. The phantom was made with 300 mL of milk at room temperature (2% reduced, Giant Eagle, Pittsburgh, PA), 400 mL of water, and 25 μL of Higgins black India ink (Chartpak, Inc., Leeds, MA). The India ink was used as the absorber and the milk was used for mimicking moving scatterers in the tissue [33]. Apart from the two-layer phantoms, no first layer added to the liquid phantom, referred to as P0, corresponds to a homogenous phantom, which in later sections will be referred to as the single-layer phantom.
Because the top layer solid phantoms were too thin for accurate measurements of optical properties with diffuse imaging techniques, thicker homogeneous phantoms were also made from the same batch of material. These top layer phantoms were cylindrical in shape, with a diameter of 8.2 cm and a height of 5.4 cm. For bottom layer, the same liquid phantom was used for measurement of optical properties as used in the two-layer phantom.
The optical properties of the both solid and liquid phantoms were measured at 690 nm and 830 nm using a frequency-domain NIRS system (OxiplexTS, ISS, Champaign, IL) and then converted to the DCS operating wavelength of 785 nm. First, ink absorption was assumed to be spectrally flat, hence the μ a at 690 nm and 830 nm from NIRS were averaged to obtain μ a at 785 nm for DCS. On the other hand, the m 0 s at 690 nm and 830 nm were fitted to a power law to extract m 0 s at 785 nm [37]. These optical properties' values (top layer: μ a = 0.11 cm −1 and m 0 s = 8.50 cm −1 , bottom layer: μ a = 0.12 cm −1 and m 0 s = 10.46 cm −1 ) were used to obtain D n during fitting of g 1 , and will be referred to as baseline optical properties in later sections. The optical properties' values for each layer were close to the range for human skull and brain tissue as reported in the literature (skull: μ a = 0.21-0.36 cm −1 and m 0 s = 11.90-7.70 cm −1 for λ = 674-956 nm, brain tissue: μ a = 0.17-0.21 cm −1 and m 0 s = 8-11.20 cm −1 for λ = 674-956 nm) [38,39]. During measurement, a solid phantom was placed on the liquid homogeneous phantom inside a glass beaker to mimic a two-layer geometry, as shown in Fig 2. The solid top layer was suspended by a wire scaffold (not shown in Fig 2) such that it was just touching the surface of the liquid layer at the bottom. The source and detector fibers were placed in contact with the top layer with the fibers being perpendicular to the surface. A black cloth was used to cover the beaker, blocking ambient light during measurement acquisition.

Influence of top layer thickness on correlation curves
The two-layer phantoms were measured with DCS, where the bottom layer liquid phantom stayed the same, but the top layer solid phantom was varied between different thicknesses, referred to as P0-P7 (P0 = no first layer, P1 = 2.01 mm, P2 = 3.28 mm, P3 = 3.87 mm, P4 = 4.65 mm, and P5 = 4.99 mm, P6 = 5.96 mm, P7 = 8.08 mm). Fig 3A shows a single trial of the measured g 1 (averaged over 100s) for all phantoms at a SD3 = 20 mm. Fig 3B shows the corresponding g 1 . In addition, photon count rates for each SD (SD1 = 10 mm, SD2 = 15 mm, SD3 = 20 mm, SD4 = 25 mm) across P0-P7 are shown in Fig 4A-4D. We can see that both g 2 and g 1 curves have distinct shape changes with various top layer thickness. This observation thus works as the basis for the feasibility of extracting the top layer thickness.

Two-layer vs. the single-layer model
For the homogenous phantom P0, Eq (14) with n = 1 was used to extract the ground-truth D B2 ; Eq (2) was used to extract the αD B . Fig 5 shows an example of the measured and fitted g 1 from a single trial of measurement for P0 at SD3 = 20 mm. The purple line and shaded region stem from the average and standard deviation of g 1 over 100 s per delay time. The green line with circle and orange line represents fitted g 1 using Eq (14) with n = 1 and Eq (2) respectively. The two homogenous fits (green and yellow lines) closely match each other, with the residual (measured g 1 minus fitted g 1 ) mostly centered around zero.
The measured g 1 curves from two-layer phantoms P1 to P7 were fitted to the two-layer analytical model to extract d 1 and bottom layer flow (D B2 ), and the conventional single-layer analytical model (Eq (2)) to extract single-layer flow (αD B ). Fig 6A-6D shows examples of the measured and fitted g 1 for the two-layer phantoms P1 = 2.01 mm, P3 = 3.87 mm, P5 = 4.99 mm, and P7 = 8.08 mm at SD3 = 20 mm. The purple solid lines with shaded regions represent the mean and standard deviation of the measured g 1 ; the green lines with circle and orange lines represent the fitted g 1 using two-layer and single-layer model (Eq (2)) respectively. We can see that although the fitted g 1 by the two-layer model did not describe the entire measured g 1 accurately, it still performed better than the one-layer model by its ability to change the slope with different d 1 . This is further quantified in terms of residuals as seen in Fig 6, where the two-layer fit shows smaller residuals at longer delay times compared to the single layer fit.

Simultaneous fitting for thickness and flow using a two-layer model
Two parameters-top layer thickness (d 1 ) and bottom layer Brownian diffusion coefficient (D B2 )-were extracted by fitting the two-layer model to the measured g 1 , while only the Brownian diffusion coefficient times α (αD B ) was extracted by the conventional single-layer model (Eq (2)). The S1 Table summarizes the results of these fitted d 1 , D B2 and αD B value for P0-P7 phantoms, where the average and standard deviation were derived from the ten repeated measurements. For visualization, Fig 7A shows correlations between the true and fitted d 1 at each SD distances for P1-P7. Linear fitting (MATLAB "fitlm") between fitted and true d 1 for all P1-P7 did not show a significant correlation. However, we found a significant positive linear correlation for all SD distances when fitting up to P5 = 4.99 mm (SD1 = 10 mm: p = 0.043 < 0.050, SD2 = 15 mm: p = 0.033 < 0.050, SD3 = 20 mm: p = 0.028 < 0.050, SD4 = 25 mm: p = 0.022 < 0.050). We observed that the fitted d 1 values at larger SD distances were overestimated more, but they follow a linear trend when d 1 < 5 mm. This could indicate that for a set d 1 , the fitting accuracy depends on the SD distance and top layer thickness.
Since the two-layer fitting yielded D B2 and the single-layer fit αD B , a direct comparison between the values was not possible. To compare the two models, we thus obtained the ground truth D B2 and αD B from measurements on P0 (homogeneous liquid phantom that provides the most accurate prediction of flow), and then calculated the percent differences between values extracted from each layered phantom (P1-P7) to their ground truths. For P0, the ground truth D B2 and αD B were calculated using single-layer models Eq (14) with n = 1 and Eq (2) respectively.
As shown in Fig 7B, for the single-layer model (triangles), fitted αD B at all SD distances quickly deviates from the baseline flow of P0 as d 1 increases (around 65% difference at the thinnest P1 = 2.01 mm), confirming its inaccuracy for flow estimation at non-negligible top layer thicknesses. On the other hand, from the results of the two-layer model (circles), we observed that the estimation of D B2 changes with SD distances: At shorter SD distances (SD1 = 10 mm and SD2 = 15 mm) fitted D B2 decreases in a similar way to the single-layer model but with smaller discrepancies compared to the reference. At longer SD distances, however, fitted D B2 only decreases to less than or around 15% from baseline until d 1 = 4.65 mm for SD3 = 20 mm and d 1 = 4.99 mm for SD4 = 25 mm. For even the thickness P7 = 8.08 mm, D B2 estimated from SD4 decrease to around 65% from baseline, which is a similar discrepancy to single-layer model at the thinnest P0 = 2.01 mm. Longer SD distances (yellow and green circles in Fig 7B for SD3 and SD4 respectively), which are more sensitive to the bottom layer, did provide better fit for D B2 with the two-layer model. These findings suggest that the two-layer model can estimate the bottom layer flow more accurately than the single-layer model when a top static layer is present.
To further demonstrate the feasibility of simultaneously extracting the top layer thickness and flow and the fitting quality, the contour plots of RSS (residual sum of squares) between the measured and fitted g 1 were generated by varying the values of d 1 and D B2 pair. Fig 8A-8D shows the example RSS for P1 (2.01 mm), P3 (3.28 mm), P5 (4.99 mm), and P7 (8.08 mm) at SD3 (20 mm). We can see that RSS in all cases converged to a single minimum as shown by the red cross. As discussed in the previous session, the accuracy of fitted d 1 and D B2 depend on

Model sensitivity on layer optical properties
The two-layer model assumes optical properties of the top and bottom layers during fitting of d 1 and D B2 . Precise knowledge of the optical properties in each layer, however, can be hard to obtain for some applications [40]. To evaluate the influence of layer optical properties on the reconstruction of d 1 and D B2 , we separately varied the μ a and m 0 s in each layer by ± 20% and performed the same two-parameter fitting procedure as described previously. The fitted d 1 and D B2 were then compared to their expected values by 100%�(d 1,fitted −d 1,true )/d 1,true and 100%�(D B2,fitted −D B2,true )/D B2,true respectively. The ground truth D B2 were fitted by the baseline μ a and m 0 s using the Eq (14) with n = 1 from P0 (liquid layer only, no top layer). For d 1 , Fig 9A-9D and 9E-9H show how changes in μ a and m 0 s affect the fitted d 1 respectively at each top layer thickness from SD1 to SD4. The overlapping curves in the figure indicate that changes in μ a and m 0 s by ± 20% in each layer minimally affect the fitted d 1 . Fig 10A-10D shows how changes in μ a affect the fitted D B2 at each top layer thickness from SD1 to SD4. We observed that changes in the first layer μ a by ± 20% (red and blue) only minimally affect the fitted D B2 . Changes in the second layer μ a , however, show that overestimation of μ a (yellow) increases the fitted D B2 by around 15%; underestimation of μ a (green) decreases the fitted D B2 by around 15%. Fig 10E-10H shows how changes in m 0 s by ± 20% affect the fitted D B2 at each top layer thickness from SD1 to SD4. We can see that for m 0 s in both layers, overestimation of m 0 s (blue and yellow) decreases the fitted D B2 ; underestimation of m 0 s (red and green) increases the fitted D B2 . In addition, changes in the second layer m 0 s affect fitted D B2 more than the first layer m 0 s (around 20% compared to around 10%).
Using the two-layer model, we found that the uncertainties in optical properties affect fitted D B2 more than fitted d 1 . Uncertainties in optical properties of ± 20% minimally affected the fitted d 1 , but those in each layer affected the fitted D B2 differently as discussed above.

Discussion and conclusion
In this manuscript, a two-layer analytical model for CW-DCS adapted from Li et al. was adapted and validated for its feasibility in estimating layer thickness in addition to flow [22]. Two-layer phantoms were fabricated to mimic the brain geometry-with the top layer representing skull with finite thickness (up to~8 mm) and the bottom layer representing brain cortex with semi-infinite geometry. DCS measurements were performed on phantoms with different thicknesses (from P0 = 0 mm to P7 = 8.08 mm). We found that, the shapes of g 2 and g 1 change with d 1 , indicating the possibility of extracting layer thickness beside flow from the DCS measurements. This finding was in compliance with Boas et al. who observed, while simulating burned tissue, a slower rate of decay in g 1 as the thickness of the teflon layer increased resting above an intralipid liquid layer [11].
The two-layer model was then used to fit for d 1 and D B2 simultaneously from the g 1 curves. For a fixed two-layer phantom, the accuracy of fitted d 1 and D B2 depends on the SD distance and the top layer thickness. For example, as shown in Fig 7A, fitted and true d 1 showed significant positive linear relationship up to P5 = 4.99 mm. In addition, when the SD is not sensitive enough to the bottom layer, as shown in the case of SD1 = 10 mm for P7 = 8.08 mm, fitted d 1 will largely deviate from its expected value. In another word, P7 = 8.08 mm is almost equivalent to a semi-infinite phantom to SD1 = 10 mm, causing inaccuracy in fitted d 1 from the twolayer model.
For the flow in the bottom layer (αD B and D B2 ), we first calculated their expected values by fitting the single-layer models (Eq (2) and Eq (14) with n = 1) to the measurements from the homogenous phantom P0 (no top layer). Then, as shown in Fig 7A for the conventional single-layer model, the percent difference between the fitted αD B quickly deviates from its expected value when solid layers were added. On the other hand, as shown in Fig 7B, fitted D B2 from the two-layer model are more stable. Especially for the longer SD that are more sensitive to the bottom layer. The longest SD4 = 25 mm produced consistently accurate fitted D B2 (less than or around 10% difference from the expected value) for P1 = 2.01 mm to P5 = 4.99 mm. D B2 deviated more from the baseline at P6 = 5.96 mm and P7 = 8.08 mm, which we attribute partially due to low photon count rate as shown in Fig 4D. To test the quality of two-parameter fitting, similar to Li et al. and Dong et al., contour plots of RSS between the fitted and measured g 1 (Fig 8A-8D) were generated as a function of the fitted parameters d 1 and D B2 ; the convergence to a single minimum demonstrated the feasibility of fitting d 1 and D B2 simultaneously [22,41].
Uncertainties in μ a and m 0 s were shown by Irwin et al. to have influence on the fitted blood flow using the conventional single-layer model [42]. To investigate how these uncertainties would influence the two-layer model, we varied the optical properties in each layer during the fitting by ± 20%. Our results (Figs 9 and 10) show that the fitted d 1 was only mildly affected within ± 20% uncertainties. In addition, μ a in the first layer only mildly affected the fitted D B2 , but overestimated and underestimated μ a in the second layer caused increased and decreased fitted D B2 respectively. On the other hand, overestimated and underestimated m 0 s in both layers decreased and increased fitted D B2 respectively. This finding is consistent with that from Li et al., Zhao et al., and Irwin et al. [22,27,42].
It is to note that, there are a few limitations of the proposed model at the current setting. For example, although d 1 and D B2 were fitted parameters, they were constrained to be within an expected range during fitting. If the first layer thickness is not within the expected range, the fitting results might be influenced. Further testing with even larger first layer thicknesses would shed light on this. Also, the accuracy of fitted d 1 and D B2 was shown to be dependent on the SD distance or the sensitivity to each layer. Thus, investigation of incorporating multiple SD during the fitting could make fitting across a variety of layer thicknesses more robust. In addition, uncertainties in layer optical properties add inaccuracy to the fitted parameters. This could be overcome by measuring optical properties with frequency domain or time domain near-infrared spectroscopy [43][44][45].
Most importantly, the two-layer model proposed in this work is undoubtedly a simplified version of human head and neglected any flow dynamics or effects of optical properties in the scalp layer. This could potentially be overcome by implementing a three-layered model as shown by Wu et al. and Zhao et al. [25,27]. However, a negligible scalp flow could also be achieved experimentally, as demonstrated by Baker et al., where an inflated cuff was placed around the head of the subject to eliminate scalp blood flow [46,47]. In such a scenario, the model proposed in this work would be applicable due to the static first layer consisting of scalp and the skull.
In conclusion, over the past decade, CW-DCS has been established as a promising new technique for continuous brain monitoring and tracking brain perfusion changes [2,13]. However, precise knowledge of extracerebral layers, such as skull, is also necessary and the reasons are twofold. First, if left uncorrected, the contribution from the skull can contaminate and underestimate cortical blood flow [28]. Second, the knowledge of continuous skull thickness could benefit detecting and monitoring skull deformation, something that is often found in patients with TBI [29]. To bridge the gap, this study adapted and validated a two-layer analytical model for CW-DCS that has the potential to sense skull thickness in addition to blood flow in the cortex layer in a continuous manner. Further studies will help improve model accuracy and widen the scope of the model for clinical applications. Upon successful implementation, this method could provide a robust and cost-effective tool for noninvasive quantification of skull thickness along with blood flow at the bedside for continuous monitoring of brain function at the clinic.
Supporting information S1